The non-linear Glasma 



Jiirgen Berges 1,2 , Soren Schlichting 1,3 

1 Institut fur Theoretische Physik, Universitdt Heidelberg 
Philosophenweg 16, 69120 Heidelberg 

2 ExtreMe Matter Institute EM MI, 
GSI Helmholtzzentrum 
Planckstr. 1, 64291 Darmstadt 

3 Theoriezentrum, Institut fur Kernphysik 
Technische Universitdt Darmstadt 
Schlossgartenstr. 9, 64289 Darmstadt 

We study the evolution of quantum fluctuations in the Glasma created immediately after the 
collision of heavy nuclei. It is shown how the presence of instabilities leads to an enhancement 
of non-linear interactions among initially small fluctuations. The non-linear dynamics leads to 
an enhanced growth of fluctuations in a large momentum region exceeding by far the originally 
unstable band. We investigate the dependence on the coupling constant at weak coupling using 
classical statistical lattice simulations for 5(7(2) gauge theory and show how these non-linearities 
can be analytically understood within the framework of two-particle irreducible (2PI) effective action 
techniques. The dependence on the coupling constant is only logarithmic in accordance with analytic 
expectations. Concerning the isotropization of bulk quantities, our results indicate that the system 
exhibits an order-one anisotropy on parametrically large time scales. Despite this fact, we find that 
gauge invariant pressure correlation functions seem to exhibit a power law behavior characteristic 
for wave turbulence. 



I. INTRODUCTION 



The great interest in relativistic heavy ion collision 
experiments is to a large part driven by its possibility 
to explore the properties of deconfined strongly inter- 
acting matter described by quantum chromodynamics 
(QCD). The past decades have revealed remarkable 
properties of the quark-gluon plasma, probably most 
strikingly its behavior similar to an ideal fluid [IJ [2]. 
However these properties are not directly accessible 
experimentally as they are encoded in the final particle 
spectra measured by the detectors at RHfC and the 
LHC [3HZ]- Consequently the extraction of medium 
properties crucially depends on theoretical input, such 
as the time when the plasma thermalizes locally, which 
has to be calculated within an ab-initio approach. 
The non-equilibrium dynamics of high energy nuclear 
collision poses a challenging problem in the underlying 
theory of QCD, which in practice can only be addressed 
with suitable approximations. In this context a field 
theoretical framework known as the 'color glass con- 
densate' has been developed, which provides a real 
time ab-initio description of nuclear collisions at high 
energies [8tiT0] . While several properties of the initial 
state right after the collision can be explored within 
this approach [TTHllj] . present studies have not yet been 
able to explain the thermalization mechanism |161 117j . 
Here quantum fluctuations may play an important role 
as they break the longitudinal boost invariance of the 
system and can be strongly amplified in the presence of 
plasma instabilities |16H33j . 

In this paper we investigate the impact of quantum 



fluctuations in the color glass condensate description 
of high energy heavy ion collisions. We work at weak 
coupling, where the presence of plasma instabilities has 
been established in previous works [T6Ul8j . and present 
results from classical-statistical lattice simulations along 
with analytic estimates. We investigate in detail the 
different dynamical stages of the system undergoing an 
instability and find that in addition to the 'primary' 
Weibel type instability [ToT - fT5] . 'secondary' instabilities 
emerge due to non-linear interactions of unstable modes. 
This mechanism is very similar to previous observation 
in non-expanding gauge-theories [El [20] and cosmo- 
logical models [M] and can be naturally understood in 
the framework of two-particle irreducible (2PI) effective 
action techniques [35] . 

This paper is organized as follows: In Sec. |TT] we 
present a short review of the dynamics of nuclear col- 
lisions in the color glass condensate (CGC) framework 
[HJ- We comment in particular on recent developments 
to include quantum fluctuations within an ab-initio 
approach [5J [TU] and show how the discussion in the 
literature is related to a classical-statistical treatment. 
In Sec. |IV| we present results from classical-statistical 
lattice simulations. We focus on non-linear effects and 
obtain the relevant growth rates and set-in times of 
primary and secondary instabilities. We find that our 
results are rather insensitive to the value of the strong 
coupling constant a s as long as a s -C I. We also 
investigate the impact of instabilities on bulk properties 
of the system such as the ratio of longitudinal pressure 
to energy density. Here our results indicate that the 
system remains anisotropic on parametrically large time 
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scales. We summarize our results and conclude with 
Sec. ED 



II. NON-EQUILIBRIUM DYNAMICS OF 
NUCLEAR COLLISIONS 

A. High-energy limit and CGC 

From a non-equilibrium point of view an ab-initio ap- 
proach to heavy-ion collisions requires to determine the 
initial density matrix consisting of two incoming nuclei 
in the vacuum and subsequently solving the initial value 
problem in quantum chromodynamics. Though this is 
beyond the scope of present theoretical methods, one 
may apply suitable approximations in the high energy 
and weak coupling limit, which make the problem com- 
putationally feasible. This is usually discussed in terms 
of the light-cone coordinates x^- = (t ± z)j ^/2, where at 
sufficiently high energies the incoming nuclei travel close 
to the light-cone, which is given by x ± = 0. The col- 
lision takes place around the time when x + = x~ = 0, 
where the center of mass of the nuclei coincides and an 
approximately boost invariant plasma is formed after the 
collision. The plasma dynamics in the forward light-cone 
{x ± > 0) is usually discussed in terms of the co-moving 
coordinates 



= v 7 ^ 



r\ — atanh(z/£) 



(1) 



where r is the proper time in the longitudi- 
nal direction and rj is the longitudinal rapidity. 
The metric in these coordinates takes the form 
9ij,v( x ) = diag(l,— 1, — 1,— t 2 ) and we denote the metric 
determinant as g{x) — det g^ v {x). The dynamics of 
the collision and the geometry of the coordinates is 
illustrated in Fig. [I] The different colors in the forward 
light-cone illustrate the dynamics of the longitudinally 
expanding plasma, which we study in this paper. 

In the color-glass framework one considers the 
dynamics of the plasma at mid-rapidity (77 <C rjBeam) 
and the nuclear partons at high rapidities separately. 
In the eikonal approximation the trajectories of the 
incoming nuclei are unaffected by the collision, while the 
dynamics of gluons at mid-rapdity is described by the 
classical action 



S[A] 



?%(x) g^{x)g^{x) K p (x) , (2) 



where J x — J d 4 Xy/—g(x) and T^ v (x) denotes the non- 
abelian field strength tensor 

T%{x) = 8Mx) - d„AZ(x)+gr bc Al(x)At(x) .(3) 

with Lorentz indices /1, v and color indices a = 1...N 2 — 1 
for SU(N C ) gauge theories. In addition, the gauge field 
A*(x) is coupled to an eikonal current Ja(x), which is 
determined by the properties of the nuclear wavefunc- 
tion at high rapidities. In practice, the separation in fast 



x =0 q = const 




FIG. 1. (color online) Cartoon of the real time evolution of 
a high energy heavy-ion collision. 



and slow degrees of freedom is performed by a renor- 
malization group procedure prescribed by the JIMWLK 
equations j8] . In the high energy limit the eikonal current 
J£(x) is given by static color sources on the light-cone 
and takes the form 

J£(t,x x ,z) = 5>»q£\x ± )8{x-) + 6»-qW(x x )5(x + ) , 

(4) 

where 5^ is the Kronecker delta in light-cone coordi- 
nates and we denote transverse coordinates as x x = 
(x 1 ^ 2 ). The color charge densities Qu^ 2 \xx), where 
the superscript (1/2) labels the different nuclei, contain 
all further information about the beam energy, nuclear 
species and impact parameter dependence. At high col- 
lider energies, these have been conjectured to exhibit a 
universal behavior, which is described by the saturation 
scale Q s and the value of the strong coupling constant 
a s [36l |37] . In this work we simply adapt the McLerran- 
Venugopalan (MV) saturation model [38] , where the color 
charge densities of the nuclei are given by uncorrelated 
Gaussian configurations 



(y±)) 



5 V 5 AB 6 ab 6(x ± - y ± ) . (5) 



The model parameter g 2 /i is proportional to the physi- 
cal saturation scale Q s up to logarithmic corrections and 
reflects the properties of the saturated wavefunctions of 
large nuclei [38]. Consequently the current in Eq. Q 
is parametrically large, i.e. formally 0(1/ g) in powers 
of the coupling constant. This makes the problem in- 
herently non-perturbative and we will come back to this 
aspect when we discuss the impact of quantum fluctua- 
tions. In addition to Eq. we impose a color neutral- 
ity constraint on the color charge densities such that the 
global color charge vanishes separately for each nucleus, 
i.e. 



d 2 x x g ( a 1/2) (x x ) = V 



(6) 
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By specifying the eikonal current according to Eq. @, 
the longitudinal geometry of the collision has effectively 
been reduced to the collision of two-dimensional sheets 
and there is no longer a longitudinal scale inherent to the 
problem. However quantum fluctuations explicitly break 
the longitudinal boost invariance of the system and may 
therefore play an important role in the non-equilibrium 
dynamics right after the collision [TO]. Before we turn 
to a detailed discussion of quantum fluctuations, we will 
briefly review the classical solution to the particle pro- 
duction process. We will show later, in Sec. |II C| how 
this solution emerges also in the weak-coupling limit of 
the quantum field theory, where quantum fluctuations 
can be handled properly. 



B. Classical solution 

Neglecting quantum fluctuations for the moment, the 
absence of a longitudinal scale in Eq. Q leads to boost- 
invariant solutions of the classical Yang-Mills field equa- 
tions 



SS[A] 

JMx) 



(7) 



In the classical color glass picture, the strong color-fields 
right after the collision are entirely determined by the 
continuity conditions on the light-cone (x^ = 0) [3"9"H4"2"] . 
Adapting the Fock-Schwinger gauge condition (A T = 0), 
where the classical Yang-Mills action takes the form 



S[A] = J 



rdr dr\ d x± 



a\2 



1 



(d T A 



a\2 



x*a -77a _ _ -77a 77a 
2 r 2 Vi Vi 4 *J ij 



(8) 



(i = 1,2) the initial state right after the collision can 
be specified at r = + , where the chromo magnetic and 
electric fields are given by 39442] 



A v =0, (9) 



Ai( Xl _) 

Ei = 0, E v (x±) = ig[aW(x±) ,4 2 \x x )] . 

Here a\ (x±) are pure gauge configurations which de- 
scribe the Yang-Mills field outside the light-cone. They 
are related to the nuclear color charge densities by [39T 



this classical evolution and leads to an effectively 2+1 di- 
mensional Yang Mills theory coupled to an adjoint scalar 
field [TlTfl"4"] . In order to study the full 3+1 dimensional 
Yang Mills theory it is therefore crucial to include quan- 
tum fluctuations, which break the boost-invariance of the 
system explicitly. 



C. Quantum fluctuations 

The inclusion of quantum fluctuations has recently at- 
tracted great attention due to their expected importance 
in understanding the thermalization process [51 HH] ■ This 
concerns in particular the inclusion of vacuum fluctua- 
tions in the initial state, which are quantum in origin 
but evolve classically at sufficiently weak coupling and 
short enough times [15] . In contrast to most discussions 
in the literature [H [TO], our analysis is based on the 
two particle irreducible (2PI) effective action framework 
[35] , This formalism has been successfully applied to a 
variety of similar problems in scalar theories [341 144j and 
gauge theories [TU [201 ES] and recently been developed 
for the problem under consideration here [46 . We give 
a short general introduction to the formalism and show 
first how the realization proposed in Ref. (TUj emerges in 
this framework. Then we go beyond that discussion and 
identify sub-leading quantum corrections and describe 
the non-linear dynamics of instabilities analytically. 

The quantum evolution equations can be formu- 
lated in terms of the expectation values of the gauge 
field operators (x) denoted by 

a°(x) = (A;(x)) , (11) 

and the time ordered two-point correlation function 



Gf u (x,y) = {TA*(x)Ai(y)) , (12) 

The two independent parts of the propagator G°^(a;,y) 
can be expressed in terms of the spectral and statistical 
two-point correlation functions 

Gf v {x,y) = F${x,y) - |sgn(x° - y°)pf v (x,y) (13) 

which are associated to the commutator and anti- 
commutator 



Xj_) 



(10) 



and dep end on transverse coordinates only. The relations 
Q and ( 10 ) specify the 'Glasma' initial state at r = + 
right after the collision. The time evolution in the for- 
ward light-cone can be studied numerically by solving the 
lattice analogue of the classical evolution equations and 
has been studied extensively [TlTH4j . However the lon- 
gitudinal boost invariance of the system is preserved in 



P Z{x,y) = i( A a Jx),Al(y) 



(14) 



= \ ({A;(x),A b u (y)}) - A;(x)A b v (y) .(15) 

Here expectation values are given by the trace over the 
initial vacuum density matrix in the presence of the 
eikonal currents. The initial density matrix is specified in 
the remote past (to — > ~°°); where the background field 
A^ vanishes and the statistical fluctuations take the 
vacuum form (see e.g. Ref. 46 ). In contrast, the initial 
values of the spectral function arc entirely determined by 
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FIG. 2. Vertices in non-abelian gauge theory in the presence 
of background gauge fields. 



equations for connected one and two-point correlation 
functions follow from the stationarity of the two particle 
irreducible (2PI) effective action g7J EES] 



T 2P1 [A,G}=S[A] 



l -t?[\nG-'] + l -tv[G- l [A]G] 
T 2 [A,G] + const . (17) 



the equal time commutation relations, which in temporal 
gauge (Aq = 0) read 



pfvfa v) 1*0=5,0 = o ) 
d x opt{x,y)\, = -8 ab 



5(x-y) , 



d x od y opf u (x, y)\ = 



(16) 



and are valid at all timesQ The gauge field expectation 
values in Eq. (11) correspond to the Glasma background 



fields, while the spectral and statistical two-point func- 
tions contain the quantum fluctuations. The evolution 



and form a closed set of coupled evolution equations. The 
set of equations is given by the evolution equation of the 
macroscopic field 



6S[A] 



-tr 



8G^ 



A] 



&A*Jx) 



G 



ST 2 [A,G] 



(18) 



and the evolution equations for spectral and statistical 
two point correlation functions, which can be written as 
(351 



iG^ a r[x;A}+n^(x)]pf u (x,y) 



dz T]V)M(x,z)p±(z,y), 



(19) 



iG^ a r{x;A}+Tl^(x)}F£(x,y) 



dz IL^(x,z)F^(z,y) + / dz U^(x, z)p*(z, y) . (20) 



Here 



G Q 'b"[x\ A] denotes the free inverse propagator 



denote dz — J dz° J d d z^/ —g(z) and 



iG, 



0,ab 



x-A] = ^(x) D?[A] 7 (s) gVgT D c a b [A] 
-^(x) D" C [A] 7 (x) gVgT D*[A] 



gr bc Fr(x)[A] 



with 7(2;) = y/—g(x) and we introduced the (back- 
ground) covariant derivative 



Df[A] 



ab 



gf abc A 



(22) 



The non-zero spectral and statistical parts of the self- 
energy I[( p / F ) [A, G] on the right hand side and the local 
part n(°)[G] on the left hand side make the evolution 
equations non-linear in the fluctuations. In general 
they contain contributions from the vertices depicted 
in Fig. [2j where in addition to the classical three gluon 
vertex there is a three gluon vertex associated with 
the presence of a non-vanishing background field. The 
explicit expressions for the derivatives on the right 



1 Note that Eq. 
placing x° = t and x A = r] and imposin, 
(A T = 0) gauge condition 



is valid also for (t, rf) coordinates, when re- 
3 — „ imt^oino- t ne Fock-Schwinger 



hand side of Eq. ( 18 ) and the self-energy contributions 
entering Eqns. (191 and (20) have been calculated to 
three loop order (g") in Ret. |45j and the corresponding 
expressions in co-moving (r, 77) coordinates can be found 
in Ref. 146 . Before we turn to a more detailed discussion 



of the right hand side contributions of Eqns. (18), (19) 



(21) and (20), it is insightful to consider first the leading part 



in a weak coupling expansion. We will see shortly how 
this recovers the classical solution for the background 
field, as discussed in Sec. |IIB| while initial state vac- 
uum fluctuations are already included to leading order 
in terms of the connected two-point correlation functions. 

In order to isolate the leading contributions one 
has to take into account the strong external currents 
J£ ~ 0(l/<7), which induce non-perturbatively large 
background fields A^(x) ~ 0(1/ g). In contrast, the 
statistical fluctuations F^ b (x,y) originate from initial 
state vacuum and are therefore initally 0(1). The 
spectral function p ab v (x,y) has to comply with the 
equal time commutation relations ( 16 1 and is therefore 



parametrically 0(1) at any time. Considering only the 
leading contributions in a weak coupling expansion, the 
evolution equation ( 18 ) reduces to its classical form (c.f. 
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Eq. @) 



SS[A] 



(23) 



and the evolution equations for the spectral and statisti- 
cal two-point correlation functions at leading order read 



iG^ a >r[x;A}p*(x,y)=0 
iG^ a >r[x;A]F%(x,y)=0 



(24) 
(25) 



where sub-leading contributions are suppressed by at 
least a factor of g 2 relative to the leading contribution. It 
is important to realize that at this order the evolution of 
the Glasma background fields decouples from that of the 
fluctuations, i.e. there is no back-reaction from the fluc- 
tuations on the background fields. Therefore the dynam- 
ics of the background fields remains unchanged and one 
recovers the classical field solutions discussed in Sec. Ill Bl 
In addition the evolution of vacuum fluctuations of the 



initial state is taken into account by Eqns. (24) and (25) 



to linear order in the fluctuations. To this order the quan- 
tum field theory is known to agree with the classical sta- 
tistical theory [13] and Eqns. (24 25 1 can equivalently be 



obtained by considering the linearized classical evolution 
equations for small fluctuations 2 |(see e.g. Ref. [10]). The 



linear approximation in Eqns. ( 241) and (25) yields a ma- 



jor simplification to Eqns. (19) and (20), as one can solve 



Eqns. (24) and (25) independently from the evolution of 
the background field. This has been exploited in Ref. [10] 
to obtain the spectrum of initial fluctuations right after 
the collision. In turn, the range of validity of the ap- 
proximation is limited to the domain where fluctuations 
remain parametrically small. This is, however, not the 
case in the forward light-cone (r > 0), where Eqns. (24) 



and ( 25 ) exhibit plasma instabilities associated to expo- 
nential growth of statistical fluctuations [HI [17] 



F^{t,t' ,x T ,y T ,v) oc exp[r(^)(^V+ y/g 2 ^)] , 

(26) 

for characteristic modes, where v is the Fourier coefficient 
with respect to relative rapidity according to 



dv 
2^ 



F%{x,y,vy^ 



-Vy) 



(27) 



In this regime fluctuations can become parametrically 
large and strongly modify the naive power counting. 
Since the approximation underlying Eqns. (24) and (25) 



is not energy conserving, one encounters exponential di- 
vergences when stressing Eqns. (24 25 1 beyond their 
range of validity [49] . In this regime it is crucial to include 



2 This can be seen immediately by solving Eq. i |2 r )[ ) in terms 
of mode functions F(x,y) = SA(x)8A(y), which then 



m- 



divdually satisfy the linearized classical evolution equations 
iG^ l [x;A] SA(x) = . 



Dynamical power counting 



W F ': one loop 



p = l 



two loop 




g 4 F 3 (p), = 4/3 



n (0) (local) 




g 2 F , = 2 



TABLE I. Self energy diagrams to two loop order (g 4 ). The 
value P corresponds to the classification in the dynamical 
power-counting scheme. The one-loop diagram in the top 
panel yields the first relevant correction, while diagrams with 
higher values of P become important at later times. 



higher order self-energy corrections, which naturally cure 
the divergencies. The associated power counting is dis- 
cussed in Sec. |IID| where we identify the diagrammatic 
contributions, which contain the first relevant corrections 
and analyze their impact on the dynamics of the instabil- 
ity. A systematic way to include self-energy corrections 
to all orders was outlined in Ref. [10] and will be discussed 
in more detail in Sec. |III| on the classical-statistical ap- 
proximation. 



D. Dynamical power counting and non-linear 
dynamics 

In order to go beyond a fixed-order coupling expansion 
of the 2PI effective action, one has to develop a power 
counting scheme which takes into account not only the 
suppression by the small coupling constant but also the 
enhancement due to parametrically large fluctuations in 
the presence of non-equilibrium instabilities. This has 
been worked out in detail for scalar theories [Ml El] and 
applies in a similar way to gauge theories [HI [5D]. In 
this power counting scheme self-energy corrections are 
classified according to powers of the coupling constant 
g as well as powers of the background field A®(x) and 
the statistical fluctuations F^(x,y). For a generic 
self-energy contribution containing powers g n F m A l p k , 
the integers n, m and I yield the suppression factor from 
the coupling constant (n) as well as the enhancement 
due to a parametrically large background field (I) and 
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parametrically large fluctuations (to). The 'weight' of 
the spectral function (fc) remains parametrically of order 
one at all times as encoded in the equal-time commu- 
tation relations. For parametrically large macroscopic 
fields A"(x) ~ 1/g one expects a sizable self-energy 
correction once fluctuations have grown as large as 
F^(x,y) ~ i/g( n - l )/ m for characteristic modes. The 
hierarchy emerging from a classification of diagrams in 
terms of ft = (n — I) jra is shown in Tab. [i] The one loop 
diagram shown in the upper panel contains two three 
gluon vertices, which give rise to a suppression factor g 2 
((n — l) = 2). On the other hand the diagram is enhanced 
by two statistical propagators in the loop (to = 2) and 
we can classify the overall contribution as 0(g 2 F 2 ). 
Similarly one can analyze the two loop diagrams and 
the tadpole diagram also shown in Tab.Jlj The tadpole 
diagram contains a suppression factor g from the four 
gluon vertex (n = 2, 1 = 0) and one statistical propaga- 
tor (to = 1), such that the overall contribution can be 
classified as 0(g 2 F). The two loop diagrams are of order 
g 4 (n — I = 4) in the coupling constant and contain at 
most three statistical propagators (to = 3). The overall 
contribution is thus 0(g 4 F 3 ) in the dynamical power 
counting. The classification in terms of /3 = (n — l)/m 
shows that the only diagram yielding a contribution to 
P = 1 is the one-loop diagram in the upper panel of 
Tab. |TJ The leading contribution of higher order loop 
diagrams can be classified as f3 = 2L/(L + 1), where 
L > 1 is the number of loops. The tadpole diagram 
yields f3 = 2. Finally there are self-energy contributions 
which contain powers of the spectral function instead 
of statistical propagators. These give rise to even 
higher values of j3 > 2 and are therefore suppressed at 
sufficiently weak coupling. 



This hierarchy of diagrams has important conse- 
quences for the dynamics of systems undergoing an 
instability, where initially small fluctuations grow expo- 
nentially in time. At early times statistical fluctuations 
are small and the system is accurately described by 
the set of equations ( 23 24j 25| , which give rise to the 
linear instability regime |16l I17j . At later times, when 
statistical fluctuations have grown larger, self-energy 
corrections become important and alter the dynamics 
of the system. In this regime self-energy contributions 
with smaller values of /3 become important at an earlier 
stage as compared to contributions with higher values 
of j3, as each diagram requires F^(x,y) ~ g~@ to yield 
a significant contribution. The one-loop diagram in the 
top panel of Tab. [i] is of order 0{g 2 F 2 ) in the dynamical 
power counting. The contribution of this diagram 
becomes of 0(1) as soon as statistical fluctuations 
have grown as large as 0(l/g), while all higher order 
self-energy corrections are still suppressed by at least a 
fractional power of the coupling constant. Accordingly 
there is a regime where the one loop diagram displayed in 
the top panel of Tab.|T]yields the only relevant correction. 

In order to investigate the impact of the one-loop 
correction in more detail, we will in the following 
neglect all contributions to the memory integrals in 
Eqns. (19 20 ) which originate from outside the for- 
ward light cone. This assumption is justified at weak 
coupling, where fluctuations are sufficiently small until 
the time when they exhibit exponential growth due to 
the presence of the instability. We can then switch to 
a discussion in (r, rf) coordinates and employ the Fock- 
Schwinger (A T = 0) gauge condition in the following. 
The self-energy contributions from the one-loop diagram 
in Tab. [I] (top panel) then take the form [IB] 



n<3#(w) = -£ 
ij£#(w) = -? 



J ^Zlc [f% V, t/)F% (x,y,v- t/) - \p% (x, y, u')p^, (x, y, «/-,/)] t*£X , (28) 
/ \F% (*, V, "')p?y (x,y,v-v') + pZ (x, y, v')F$ (x, y, v-t/)] Sf^, d , (29) 

J v' 



where x = (t x , x±) collectively labels the proper time and 
transverse coordinates, v is the Fourier coefficient with 
respect to relative rapidity and Lorentz indices take the 
values p = 1,2, -q. The three gluon vertices in Eqns. (28) 
and ( 29 1 take the form 



V: 



v: 



(30) 



"7 _ y -^vj 

x,abc x,abc v x,abc 

where in addition to the classical three gluon vertex there 
is a contribution from the background field. The classical 
three gluon vertex can be written as 



V 



x,abc 



f 



a be 



g" v (x)(-d2' : 



28^) 



(31) 



g^(x)(dr - d^) + g^(x)(2d^ + 8^ 



(no summation over b and c), where the derivative oper- 
ator d%' x = (d Tx , —d x± , ~iT~ 2 v) only acts on the propa- 
gator with color index a, and the contribution from the 
background field is given by 

VxS? = ~9 [ <T (z)C afe)Cd A}(x) (32) 
+ g^Cb cM A» d (x) + g^C caM A» d {x)] 

with 

C ab , cd = f ace f bde + f ade f bce . (33) 

In order to investigate the impact of the above self-energy 
corrections, we first note that the pp term in Eq. (28) is a 



7 



genuine quantum correction, which is absent in the clas- 
sical statistical theory [43] . On the other hand statistical 



fluctuations grow exponentially according to Eq. ( 26 ) for 
unstable modes, such that the dominant contribution in 



Eq. ( 28 ) arises from the classical (FF) contributions and 
one can safely neglect the sub-leading quantum correc- 
tions. The right hand side of the evolution equations ( 20 ) 



then receive exponentially enhanced contributions from 



the self-energies ( 28 1 and ( 29 ) , which grow exponentially 
in (proper) time. This behavior can be verified explic- 
itly by performing a 'memory expansion' of Eq. (|20|, i.e. 



evaluating the memory integrals on the right hand side 
of Eq. ( 20 1 around the latest (proper) times of interest 
[34l [44] . We find that the dominant contribution origi- 
nates from the statistical self-energy in Eq. ( 28 ) , whereas 



the contribution from the spectral self-energy in Eq. ( 29 ) 
is effectively f3 = 2 in the above classification scheme and 
thus becomes important only at later times. The modi- 
fied evolution equations in this regime then take the form 

iG^r[x,v;A]F%(x,y,v) = f C^x, y, v)g^{y) , 

(34) 

where S T is the extent in time for which the memory 
integrals are evaluated. To obtain Eq. (34 1, we per- 



formed a leading order Taylor expansion of the integrand 
in Eq. ( 20 ) around t z = t v and made use of the equal time 



commutation relations in Eq. ( 16 ) to estimate the spec- 
tral function. The one-loop integral II^) , which appears 
on the right hand side of Eq. ( 34 1 , is dominated by the 
contributions from unstable modes and is proportional 
to exp [2T(j^g 2 fiT] at equal times. The contribution on 
the right hand side acts as a source term in the evolu- 
tion equation of statistical fluctuations. As is well known 
from various examples of self-interacting quantum field 
theories, this term leads to a non-linear amplification of 
instabilities, where 'secondary' instabilities with strongly 
enhanced growth rates emerge over a large range of mo- 
menta. This has been observed in scalar field-theories 
[341 14*4] as well as in non-abelian gauge theories [T5J [22] 
and is a rather generic feature of self-interacting theories 
undergoing an instability. We will show in Sec. |IV| that 
the phenomenon of non-linear amplification also emerges 
in numerical simulations of the unstable Glasma and 
plays a crucial role in understanding gauge-invariant ob- 
servables. The characteristic time scale for non-linear 
amplification to take place can be infered by comparing 
the magnitude of the non- linear contributions in Eq. ( 34 ) 



to the contributions of the background field. In the weak 
coupling limit this time scale is paramctrically given by 



1 

2Tq 



In g~ 



(35) 



where Tq is the characteristic growth rate of primary 
instabilities and we assumed Ff v ~ 0(1) at r = 0+ for 
characteristic modes. In addition to Eq. (351 there are 



sub-leading contributions associated to the delayed set 
in of primary instabilities and the spectral distribution 



of statistical fluctuations in the initial state. The prior 
give rise to a constant contribution -\/(? 2 /iTp r j mary , while 
the latter enter only logarithmically in this estimate. 

The emergence of secondary instabilities again modifies 
the power-counting, and one has to take into account 
also the contributions which originate from modes 
which exhibit secondary instabilities. Also higher order 
self-energy corrections become increasingly important as 
time proceeds and non-linear amplification can repeat 
itself, until at some point the growth of instabilities 
saturates and occupancies become as large as 0(l/g 2 ). 
In this regime every truncation at a fixed loop order 
breaks down and the problem has to be addressed in a 
fully non-perturbative way. While in scalar quantum 
field theories there are different ways to address this 
problem, involving e.g. large N resummation techniques 
[50] . the most frequently employed approach in gauge 
theories is the classical-statistical approximation. 



III. CLASSICAL-STATISTICAL 
APPROXIMATION 



In the classical statistical approximation all fluc- 
tuations evolve classically and one neglects genuine 
quantum fluctuations, such as the (pp) quantum term 
in Eq. (28). When transforming to the language of 



expectation values, as employed in the previous section, 
it can be shown that the prescription resums an infi- 
nite subset of diagrams [33], such that the dynamics 
of fluctuations is treated on equal footing with the 
background fields. The set of diagrams included in the 
classical-statistical treatment can be identified as the 
self-energy corrections, which contain the most powers of 
the statistical propagator as compared to powers of the 
spectral function for each topology [43 . Accordingly, 
this corresponds to resumming the leading effects of 
the instability to all orders in the coupling constant 
[10] . In contrast to expansions at fixed loop orders, 
the classical statistical approximation thus provides a 
robust approximation scheme, which is particularly well 
suited for problems involving large statistical fluctua- 
tions. However there are problems associated with the 
Rayleigh- Jeans divergence, which concern the handling 
of ultra-violet divergences and the approach to thermal 
equilibrium at late times, which are discussed in more 
detail in the literature [33] , 

In the classical-statistical theory observables (0(x)) cs 
are calculated as an ensemble average of classical 
field solutions A c \[A To , E To ], which individually satisfy 
the Yang-Mills evolution equations. The canonical 
field variables A To and E Ta at initial time To are dis- 
tributed according to a phase space density functional 
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W[A To , E To ], such that 



(0(x)) cs = J DA Ta DE T0 W[A T01 E Ta ] O cl [A To , E TQ ] 



Here O c \ [A Tg , E TQ ] denotes the fact that the observable is 
evaluated as a functional of the classical field solution 

O c i[A T0 ,E T0 ] = J DA 0[A] S(A-A cl [A T0 ,E T0 }) , 

(37) 

is the classical Yang-Mills field so- 



where A cX [A TQ , E To ] 
lution with initial conditions A c 



A Tn and E, 



m 



E T 



at initial time r . In practice, Eqns. (36) and (37) state 
that vacuum fluctuations of the initial state are added 
on top of the background field at initial time, while 
the subsequent classical evolution keeps track of all non- 
linearities. The set of equations (36) and (37) is precisely 
the same as in Ref. [TDJ, where it has been obtained as 
a partial resummation scheme of the perturbative cor- 
rections due to vacuum fluctuations of the initial state. 
The procedure outlined in Ref. |10j consist of a hybrid ap- 
proach, which employs the linearized evolution equations 



( 23 24 25 1 outside the forward light-cone, where fluctu- 
ations are small, while switching to a classical-statistical 
description in the forward light-cone. This has the ad- 
vantage that the evolution of the background field and 
the spectrum of fluctuations, which enter the phase-space 
weight W [A TQ , E To } , can be obtained analytically. The 
phase-space average in Eq. ( 36 1 can then be taken on the 



Cauchy surface tq = + , such that only the dynamics in 
the forward light-cone has to be studied within classical- 
statistical lattice simulations. We follow this approach 
but employ a simpler spectrum of initial fluctuations, as 
specified in Sec. |III B| instead. This simplification is jus- 
tified at sufficiently weak coupling, where the spectrum 
of fluctuations is quickly dominated by the growth of pri- 
mary instabilities rather than the initial spectrum. 



A. Coupling dependence 

An important property of the classical-statistical de- 
scription is the independence of the gauge coupling con- 
stant g, in the sense that the classical-evolution equations 
are invariant under a change of variables 



A a Jx) g A"Jx) 



Ff v (x,y)^g 2 F^(x,y),(38) 



whereas the spectral function remains unaffected, as en- 
coded in the equal time commutation relationsr] With 
the rescaling ( 38 ) the entire coupling dependence in the 



classical-statistical evolution can be absorbed into the 



3 In the classical theory —i times the commutator is to be replaced 
by the Poisson brackets 43 . 



initial conditions at r = + , while the classical evolu- 
tion equations become independent of the coupling con- 
stant. For the Glasma background fields this can be 
achieved most efficiently by replacing gA®(x) — > A®(x) 



(36) and gg^ 2 Hx ± ) 



qW 2 ){x±) for all 



expressions m 



Sec. II B where in the MV model prescription 

(~gi A Hx ± )g(B) b (y±)) = gVS{x± y±) ■ (39) 

Here the model parameter g 2 \i is directly related to the 
saturation scale Q s without further powers of the cou- 
pling constant appearing in the expression |38j . Accord- 
ingly the defining equation ( 39 1 is indeed independent of 
the value of the coupling constant. In contrast, the cou- 
pling constant g appears explicitly in the initial spectrum 
of fluctuations, given by 



=0+ 



= g 2 Ffi(x >y )\ T=T , = 



o+ 



(40) 



and similarly for derivatives at initial time. Here it is 
important to note that the magnitude of the vacuum 
fluctuations F^(x,y) on the right hand side is indepen- 
dent of the value of the coupling constant. Thus the 
initial suppression of vacuum fluctuations compared to 
the boost invariant background fields is the only mea- 
sure of the strong coupling constant, present in the 
classical-statistical field theory. We will exploit this fact 



in Sec. IV B where we vary the amplitude of initial fluc- 
tuations in our simulations to study the coupling depen- 
dence of our results. 



B. Initial conditions at r = + 

Within the classical-statistical framework, the initial 
conditions for the time-evolution in the forward light- 
cone are given at r = + by the set of equations |9j 
and (10), complemented by the spectrum of initial fluc- 
tuations. While in general it is necessary to implement 



the spectrum of fluctuations as specified in Ref. [TD], it 
is also clear that at sufficiently weak coupling many de- 
tails of the spectral shape of the initial fluctuations be- 
come irrelevant due to the presence of instabilities, which 
quickly dominate the spectrum. Since implementing the 
spectrum of fluctuations in Ref. |10| is numerically chal- 
lenging, we will therefore stick to a simpler choice, where 
in accordance with previous works [TDJ [T7] the statistical 
fluctuations initially take the form 

' , (41) 





t—t' 


=0 


d T Ff v (x,y) 


t—t' 


=0 


d T d T ,Ff„(x,y) 


t—t' 


=0 



with 



6Ef(x ± ,r,) 
6E v (x±,T)) 



d v f(v)eUx±) 
-f(ri)D i ^(x ± ) 



(42) 
(43) 



such that SE®(x) is an additive contribution to the back- 
ground field E®(x) at initial time. The advantage of this 
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FIG. 3. (color online) Time evolution of the pressure-pressure correlator TIl(t, v) for different rapidity wave numbers v. 
Once the initial fluctuations have grown larger one observes the emergence of the secondary instabilities, with growth rates. 
Subsequently the instability propagates towards higher momenta until saturation occurs and the system exhibits a much slower 
dynamics. The results are obtained for the MV model parameter g 2 fj,L± = 22.6 on latt ices with N± = f 6 and N v = f 024 sites. 
The initial fluctuations are parametrized by A = 10 -10 and b — O.Of according to Eq. (44 1. 



construction is that the Gauss constraint is satisfied ex- 
plicitly for arbitrary functions f{rj) and ef(x±). For the 
numerical simulations we choose 

{eliPT^iqT)) = A 2 S. l3 S ab 6(p T + q T ) , (44) 



{f{v)f{v'))=e-^6{v + v'), 



(45) 



where pr and v are the Fourier coefficients with respect 
to relative transverse coordinates and relative rapidity. 
Here & is a (small) number, which regulates the ultra- 
violet divergence and the dimensionless parameter A con- 
trols the initial amplitude of small wave-number fluctu- 
ations. While this construction does not respect the de- 
tails of the spectral composition, the parameter A pro- 
vides a measure of the coupling constant A 2 ~ g 2 (see 



Sec. Ill A) and we will vary its size in Sec. IV B to study 



the coupling dependence. 



IV. LATTICE RESULTS 

In this section we present results from classical- 
statistical lattice simulations of the Glasma evolution in 
the presence of boost non-invariant fluctuations. While 
the existence of a non-equilibrium instability has been 
established in previous simulations [16j E], we focus 
on the non-linear regime where unstable modes have 
grown large enough to significantly alter the dynamics. 
In contrast to the linear regime, where the initial size 



of boost non-invariant fluctuations is irrelevant for the 
dynamics of unstable modes, it is clear that for the 
non-linear regime the size of the initial fluctuations 
matters. In view of the spectrum of initial fluctuations 
obtained in Ref. [10] . the ratio of the initial amplitude of 
fluctuations compared to the amplitude of the (squared) 
background field is parametrically of the order of the 
strong coupling constant g 2 . We study this dependence 
in Sec. |IVB| by considering different amplitudes of the 
initial fluctuations, as characteri zed by the dimensionless 



parameter A ~ g (see Sec. Ill B ) . We restrict our 



analysis to weak coupling (g 2 -C 1), where classical- 
statistical methods are expected to provide an accurate 
description of the quantum dynamics on large time 
scales. 

The discussion of our results is organized as fol- 



lows: In Sec. IV A we investigate the dynamics of the 
instability at weak coupling and show how deviations 
from the linear regime emerge in terms of secondary 
instabilities. To further analyze this behavior, we obtain 
the relevant growth rates for primary and secondary 
instabilities as well as the corresponding set-in times. 
Subsequently, in Sec. |IVB[ we investigate the depen- 
dence on the coupling constant by varying the size of 
initial fluctuations. 

If not stated otherwise we perform simulations on 
N± = 16, N v = 1024 and N± = 32, N n = 128 lattices 
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and we employ the set of parameters g 2 fi N±a± = 22.6 
and N v a v = 1.6 in accordance with Ref. [16]. We study 
the time evolution of the gauge-invariant pressure- 
pressure correlation function 



n!(T,W) = (P L (T,x ± ,r))P L (T,yx,V% 



(46) 



where Pl{x) is the longitudinal pressure as a function 
of space and time arguments and < ., . >t denotes av- 
erage over transverse coordinates and classical-statistical 
ensemble average. We will frequently employ the Fourier 
transforms of H 2 l {t, r), rf) with respect to relative rapid- 
ity, i.e. we consider 

I1 2 l (t,v) = (N^y 1 J dr, dr,' n L (r,V,v') , 

(47) 

and usually show results for Hl(t, v), i.e. the square root 
of the above expression. The details of our lattice setup 
are described in more detail in the appendix. 

A. The unstable Glasma 

In a first step we study the time evolution of the 
gauge-invariant pressure-pressure correlator IT^ (t,i>) 
for a fixed value of the amplitude of initial fluctuations 
A = 10~ 10 and b = 0.01. The results are shown in Fig. [3] 
for different rapidity wave numbers v as a function of 
time. From Fig. [3] one observes a sequence of differ- 
ent dynamical regimes which are characterized as follows: 

At very early times y/g 2 ^T < 2 one observes a pe- 
riod of rapid initial growth, which is presumably caused 
by the dephasing dynamics of the strong background 
fields, that takes place roughly on the same time scale 
[TlTfl"4l [T71 1ST] . However at weak coupling, i.e. for small 
fluctuations, this constitutes a rather small effect as the 
unstable modes exhibit their dominant growth at later 
times. 

The rapid initial period is followed by a regime 
where the Glasma instability [16j [17] is operative and 
modes with non-zero rapidity wave number exhibit 
exponential amplification. The instability sets in with a 
delay for higher momentum modes and the functional 
form is well described by an exponential of the form 
exp[r(z/)-\/.9 2 A i ' r ], with the momentum dependent growth 
rate T(i/), as seen for v = 4, 12 in Fig. [3] To further 
investigate this behavior we fit a set of continuous 
piecewise linear functions to the modes displayed in 
Fig. [3] in order to obtain the relevant growth rates and 
set-in times. The results of these fits are shown in Fig. [4] 
as a function of rapidity wave number v. From the upper 
panel of Fig. [4] one observes that the primary set-in 
times follow a linear behavior, as reported in Ref. [To] . 
The primary growth rates are shown in the lower panel 
of Fig. [4] One observes that modes with small rapidity 
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FIG. 4. (color online) (top) Set in times of primary and 
secondary instabilities as a function of rapidity wave number 
v. Once secondary instabilities set in, the growth quickly 
extends to higher rapidity wave numbers, (bottom) Growth 
rates F(i>) as a function of rapidity wave number v for primary 
and secondary instabilities. The results are obtained for the 
MV model parameter g 2 /iL± — 22.6 on lattices with N± = 16 
and N„ = 1024 sites. 



wave number exhibit smaller growth rates as compared 
to modes with higher rapidity wave number, while at 
large v the primary growth rates become approximately 
constant. The numerical values are compatible with the 
results reported in Ref. [IB], where characteristic growth 
rates were obtained from a convolution of the spectrum. 

While the primary instability continues to set-in 
for higher momentum modes, one observes from Fig. [3] 
that at later times modes with intermediate (y = 43, 71) 
and small (y = 4) rapidity wave number suddenly 
exhibit much higher growth rates than previously 
observed. This change in the dynamics becomes evident 
when shortly after modes with even higher rapidity 
wave numbers (y = 94, 200) exhibit even stronger 
growth rates, such that the spectrum extends quickly to 
the ultra-violet and the instability propagates towards 
higher momenta. This is precisely the signature of 
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secondary instabilities, where non-linear self-interactions 
among unstable modes give rise to an amplification 
of the primary instability. The amplification happens 
initially in a small momentum region and then quickly 
propagates outwards to higher momenta. This can be 
seen in the upper panel of Fig. |4j where we show the 
set-in times of primary and secondary growth. One 
also observes that for modes with large rapidity wave 
number v > v c secondary instabilities set-in before 
the primary instability, such that the growth of high v 
modes is solely due to non-linear effects. The numerical 
value of v c depends, of course, on the size of the initial 
fluctuations and we will confirm the non-linear origin 
of this phenomenon in Sec. |IVB[ where we investigate 
in more detail the dependence on the initial amplitude 
of fluctuations. The dynamical power counting scheme 
developed in Sec. |IID| suggests that these secondary 
instabilities are caused by the one-loop diagram shown 
in the upper panel of Tab. [TJ with secondary growth 
rates as large as twice the primary ones. If we compare 
the rates of primary and secondary growth, as shown 
in the lower panel of Fig. |4j we find that this is indeed 
the case for modes with intermediate z/, which exhibit 
the earliest non-linear amplification. Modes with higher 
values of v exhibit even larger growth rates, which can 
be attributed to multiple amplification processes as well 
as higher order corrections. 

The growth of primary and secondary instabilities 
in Fig. [3] continues until at some point saturation of 
the instability sets in and the system has reached non- 
perturbatively large occupation numbers. In this regime 
we observe that the process of non-linear amplification 
continues even after the growth of the leading primary 
modes has saturated. This has a significant impact also 
on bulk observables such as the ratio of longitudinal 
pressure to energy density. Before we turn to a more 
detailed discussion of this highly occupied regime, we 
will first investigate the coupling dependence of the 
non-linear amplification process. 



B. Coupling dependence 



In Sec. IIV Al we have discussed the time evolution 
of non-boost invariant fluctuations in the Glasma. 
We have shown that, at weak coupling, primary in- 
stabilities of small rapidity wave number fluctuations 
occur. In turn, these cause secondary instabilities 
of modes with higher rapidity wave numbers until 
saturation of the instabilities occurs and one enters 
the highly over-occupied regime. In this section we 
investigate in more detail the dependence on the choice 
of parameters, in particular the impact of the initial 
amplitude of fluctuations. According to Sec. Ill A 



this can be interpreted as varying the value of the 
coupling constant g 2 in our simulations, without respect- 
ing in detail the spectral shape of initial fluctuations [TU] . 
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FIG. 5. (color online) Set in times of secondary instabilities 
as a function of the size of initial fluctuations. The parame- 
ter A quantifies the magnitude of initial fluctuations and is 
proportional to the coupling constant g (see also Sec. Ill B I . 
One observes a logarithmic dependence in accordance with 
analytic expectations. The results are obtained for the MV 
model parameter g 2 [iL± — 22.6. 



In order to study the dependence on the initial 
amplitude of fluctuations, we vary the parameter A in 
the range of 10 -15 to 10 -5 . The qualitative behavior is 
the same as observed in Sec. |IV A| for all values of A, i.e. 
we observe primary instabilities followed by non-linear 
amplification and subsequent saturation of the growth. 
Due to the non-linear origin, the time scales for the 
set in of secondary instabilities and the saturation of 
growth depend, of course, on the initial amplitude of 
fluctuations. The results of our analysis are shown in 
Fig. [5j where we show the characteristic set-in times 
of secondary instabilities as a function of the initial 
amplitude A of boost non-invariant fluctuations. The 
rather weak dependence observed in Fig. [5] stems from 
the fact that, at early times, the magnitude of non-linear 
contributions depends exponentially as exp[2r -\/(7 2 /ir], 
whereas the dependence on the initial amplitude is just a 
power. Assuming that non-linear amplification is caused 
by the one-loop diagram depicted in Tab. [I] we obtain 
the parametric estimate (see Sec. II D ) 
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1 

2T\ 



(48) 



where Tsetin characterizes the set-in time of primary in- 
stabilities and r is the characteristic primary growth 
rate. This behavior is reproduced by the lattice data on 
a qualitative level. 



C. Saturated regime 

The evolution of the system in the saturated regime is 
of great interest, when studying the thermalization pro- 
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FIG. 6. (color online) Ratio of longitudinal and trans- 
verse pressure as a function of time. The different curves 
correspond to different values of the lattice spacing. One ob- 
serves a remaining order one anisotropy over a large time 
scale. The results are obtained for the MV model parame- 
ter g 2 fj,L± = 22.6 and the initial fluctuations are chosen as 
A = 1(T 5 and 6 = 0.1. 



cess at weak coupling. In this regime the system exhibits 
a much slower dynamics and one expects the system to 
become approximately isotropic on sufficiently large time 
scales. We first analyze this behavior, by considering 
the evolution of the ratio of longitudinal pressure to en- 
ergy density, as a measure of the bulk anisotropy of the 
system. The challenge in this analysis comes from the 
fact that the relevant ultra-violet cutoff associated with 
longitudinal momentum A z ~ 7r/(ra^) decreases with 
(proper) time. Furthermore having a large rapdity cutoff 
~ Tr/a v can cause severe problems at early times and a 
proper renomalization scheme might be needed to ensure 
physical results. We adress this problem by chosing the 
initial amplitude of fluctuations very small, such that the 
overall contribution of fluctuations to the energy density 
is less than a percent even for the largest cutoffs that we 
consider. We then vary the lattice spacing a„ while keep- 
ing NrjCLri fixed to study the sensitivity to the cutoff. The 
results are presented in Fig. [6j where we show the ratio 
of longitudinal and transverse pressure as a function of 
time. While at early times the longitudinal pressure of 
the system is consistent with zero, we observe a clear rise 
of the longitudinal pressure towards later times. In the 
saturated regime the trend towards isotropization slows 
down dramatically and the system exhibits a remaining 
order one anisotropy over a large time scale. The results 
are insensitive of the longitudinal discretization, as long 
as the lattice spacing a v is sufficiently small. 

We also study the spectrum of the pressure-pressure cor- 
relation function II(t, v) in this regime. This is shown 
in Fig. [7J at different times, as a function of longitudi- 
nal momentum p z = v/t. From the top panel of Fig. [7J 
one observes how the spectrum rapidly extends towards 




FIG. 7. (color online) Spectrum of the pressure-pressure 
correlation function II (r, u) as a function of longitudinal mo- 
menta p z — u/t at different times of the evolution. The top 
panel shows the spectrum on a log-linear plot, whereas the 
lower panel corresponds to a double logarithmic plot. The 
black dashed line corresponds to a fit of the functional form 
(491. 



higher momenta at times g 2 /J.r ~ 1000 — 2000, while the 
amplitude at low momenta decreases. This redistribu- 
tion is accompanied by the increase in the bulk pressure 
observed in Fig. [6] At later times g 2 /iT > 2000 the evo- 
lution of the spectrum features enhanced contributions 
from soft modes and a strong fall-off at high momenta. 
The evolution of the spectrum in this regime proceeds 
in a much slower way. The lower panel of Fig. [7J shows 
the spectrum on a double logarithmic plot. One observes 
that the soft tail of the spectrum can be described by a 
power-law. This is illustrated by the black dashed line in 
Fig. [JJ which corresponds to the functional form 

f{x) = A x~ p /{l + exp[Bx 2 }) , (49) 

i.e. a power-law with normal UV regulator. This behav- 
ior is very similar to wave-turbulence observed in non- 
expanding systems [57H59j . The power law features an 
infrared power-law exponent of ft ~ 1/3, while the regu- 
lator controls the fall-off at high momenta. The shape of 
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the spectrum at later times is qualititive similar. How- 
ever we already find sizable amplitudes for modes with 
longitudinal momenta on the order of the transverse lat- 
tice cutoff. To avoid discretization errors, one therefore 
has to consider also much larger lattices in the transverse 
direction and we expect future simulations to extent the 
studies of this regime. 



V. CONCLUSION 

In this paper we discussed the impact of quantum 
fluctuations on the non-equilibrium dynamics of the 
Glasma. The picture that emerges at weak coupling is 
quite universal and depends only weakly on the details 
of initial fluctuations. In the weak coupling scenario 
initially boost non-invariant fluctuations exhibit expo- 
nential growth until at some point they have grown 
large enough for non-linear interactions to become 
important. At this stage secondary instabilities set in for 
a large momentum region, which extends to much higher 
rapidity wave numbers as the primary growth. Both 
primary and secondary growth continues until saturation 
of the instability occurs and the system exhibits a much 
slower dynamics. This scenario is similar to various 
examples of instabilities in self-interacting quantum field 
theories (THl EM EH EI] and we have shown in Sec. |IID 
how the emergence of secondary instabilities can be 
studied systematically within the framework of two 
particle irreducible effective action techniques. At the 
qualitative level this analysis is sufficient to predict the 
parametric dependence of the set-in time of secondary 
instabilities and the estimate can be made quantitative 
when the microscopic dynamics of the primary instabil- 
ity is described analytically. The set in time of secondary 
instabilities depends on the growth rates of the primary 
instability and to a much weaker extent on the size of the 
initial fluctuations, which is related to the value of the 
strong coupling constant. We confirmed this qualitative 
behavior by varying the size of the initial fluctuations 
in the classical-statistical lattice simulations, without 
respecting in detail the spectral composition. The latter 
will be taken into account in future studies, though one 
does not expect qualitative changes at weak coupling, 
where the longitudinal spectrum is quickly dominated 
by the exponential growth of primary instabilities. In 
contrast, changes in the spectrum of the background 
field may have a much more significant impact on the 
dynamics, as they readily alter the dynamics of primary 
instabilities. In particular this may change the character 
of the instability from the Weibel type, which is present 
in the MV model, to the Nielsen Olesen type jSU [52] and 
it will be important to consider more realistic saturation 
models in the future. This is important also in view 
of the applicability at RfflC and LHC energies, where 
in addition one has to consider much larger values of 
the strong coupling constant. This is not unambiguous 
since one encounters conceptual problems concerning 



the renormalization of the theory as well as the impact 
of sub-leading quantum corrections. While exploratory 
studies in scalar field theories have recently obtained 
promising results [53J [53], we expect future studies of 
non-abelian gauge theories to expand on this issue. 

In addition to the unstable regime, we also studied 
the dynamics of the saturated regime, which is of great 
interest in the recent debate on the thermalization 
mechanism at weak coupling and high collider energies 
[551 [55] , In this context, different scenarios involving 
instability induced isotropization [55) as well as the 
formation of a transient condensate |56j have been 
proposed and are currently under investigation. While 
for non-expanding systems, the results from classical- 
statistical lattice simulations suggest the occurrence of a 
turbulent cascade [57 59 , there are only few results for 
expanding systems |17[ 133) and we expect more studies 
in the future. While a dedicated lattice study has the 
potential to clarify these questions, we restricted our 
analysis to the characteristic properties of the system 
at later times. We investigated the ratio of longitudi- 
nal and transverse pressure as a measure of the bulk 
anisotropy of the system. This observable can be used 
to distinguish between different scenarios which predict 
a characteristic evolution of this quantity. Our results 
indicate that the system exhibits an order one anisotropy 
on large time scales, which appears to be in contrast to 
the scenario developed in Ref. [55], where the authors 
considered parametrically small and large anisotropies. 
By investigating the spectrum of the pressure-pressure 
correlation function, we found evidence for a scaling 
solution and similar results have also be obtained in 
Ref. [T7], where a different correlation function has 
been considered. The question whether this behavior is 
indeed related to the emergence of a turbulent cascade 
as observed in non-expanding systems [57H59_], will be 
subject to future studies. 
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APPENDIX: CLASSICAL-STATISTICAL 
LATTICE SIMULATIONS 

The classical-statistical lattice simulations are per- 
formed in the Hamiltonian framework with the lattice 
link variables U^(x) which are related to the continuum 
variables by 



U l (x) = exp[igaj_ t a A°;(x)] , 
U v {x) = exp[iga v t a A%(x)] , 



(50) 
(51) 



where a± and a n are the lattice spacings in the transverse 
and longutidinal directions and t a denote the generators 
of the SU(2) gauge group. The SU(2) exponential can 
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be computed as 

cxp[it a A a ] = cos(VI^/2)lL+ Sm( ^ /2) A a r .(52) 

v A 2 

In Fock-Schwinger gauge one finds U T = 1L and the (di- 
mensionless) electric field variables read 



E?(x)=ga rd T A?(x) , 
E«(x)=ga 2 r^drA^x) 



(53) 
(54) 



Without loss of generality we set the coupling constant 
g = 1 in the following (see Sec. Ill A). The (dimension- 
less) lattice Hamiltonian density H(t, x± , rf) of the sys- 
tem is given by 



f 2 F 2 



8| 



B\ 



(55) 



2 2 

where ce — a 2 /t 2 and c B — a 2 /(r 2 a 2 ), and the squares 
of the electric and magnetic field strengths are, for the 
SU(2) gauge-group, given by 



eUx) = [E?{x)f 



Sl(x) = [E«{x)f 



(56) 



B T (x) = 4 tr 11- U°{x) , B L {x) = 4 tr 11- U^{x) 
where U a p(x) are the standard plaquettes given by 

U^(x) = U a (x)U {x + a)Ui{x + P)Ul(x) (57) 
V^(x) = U a (x)ul(x + &- p)Ut(x - $)U p (x - $) 

where a,/3 = 1, 2, 77 and we also defined the V^L plaque- 
ttes for later use. 

A. Evolution equations 



The Hamiltonian evolution equations on the lattice in 
co-moving (r, rj) coordinates are given by the equation of 
motion for the link variables 

U a (x T ,ri,T + o T ) = W a (x T ,r),T + a T /2)U a (x T ,r),T) , 

(58) 

(no summation over a = 1,2, rf). Here denote pla- 
quettes involving a time link, given by 



W°{x) = exp [it a c a E a a {x)} 



(59) 



(no summation over a = 1, 2, r/) and the coefficients read 

(60) 



a T 

T 



Ta v a 
a 



2 



The update rules for the chromo-electric fields are given 
by 



E a a (x T , V, T + Or/2) = E a a (XT, V, T - Or/2) 



(61) 



2 d a p tr 



it a (U°(x) + V°) (a 



with the coefficients 



a- 



a T 

ra r , 



(62) 



The time evolution can then be computed by alternately 
solving Eqns. (f58l and (61 ). The Gauss law constraint is 



conserved by this evolution and reads 
a 

a=l,2,i) 



T 1 o „ 







(63) 



which is satisfied separately for all x and a. Here we 
denote hi — 1, hd v = a 2 /(r 2 a 2 ) and 



02(a:)=tr zi a W°(cc) 



(64) 



tr [rf Q [/t (a; - a)W°(a; - a)U a (x - a) 



B. Initial conditions on the lattice 

We first generate sets of uncorrelated standard Gaus- 
sian random numbers for every position in the transverse 
plane and every color associated to the color-charge den- 
sities of the nuclei, i.e. 

^V)= 9 V(^ A W-#i A) ) (65) 

where (xt) are Gaussian random number and 

(A) 

the subtraction of the overall color charge Ra = 
N±_ 2 J2 X± £a(%±) ensures the overall color neutrality 
constraint. The result is then Fourier transformed to 
momentum space where we solve the Laplace equation 



a^V) = -^ 2 p1 a V) 



(66) 



The result is Fourier transformed back to obtain the so- 
lution of the Laplace equation in coordinate space. We 
then proceed by calculating the pure gauge solutions 
J/^ 1 / 2 ' (x±) according to 



V {A) (x ± ) = exp\it a A^\x±)] . 



(67) 
(68) 



Finally the link variables U^(x±,r]) at initial time are 
obtained as [TH1 [T7] 

Ui(x x ,v) = Mi(x ± )Ni(x x ) , U n (x x ,v) = 1 , (69) 



Mi(x ± )= uP(x ± ) + Ur(xx) 



(2), 



Ni(x ± ) = 



U} {1) (x ± ) + U}M(X ± ) ^ . 
and the electric fields E"(x± ) rj) are given by [16] 
E°(x ± , V ) = tr[iP J2 U?\x ± )-U?\x ± -i) (72) 

+ uPix^ulixL) - Uj(x ± - ip( A) (x ± - 1)} . 



(70) 
(71) 
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To generate the boost non-invariant fluctuations we first 
generate the functions f{y) and ef(p±) in momentum 
space according to 



/(i/) = e -" H v^v; (73) 

e?(px) = <?a A£?(p ± ) (^Vxa) (74) 



where are uncorrelated Gaussian random num- 

bers. After performing a Fourier transform to coordinate 
space, the fluctuations 5E® (xj_, rf) are calculated as 



C. Lattice observables 

We will denote the diagonal components of the stress- 
energy tensor by the energy density e(x), the transverse 
pressure Pt{x) and the longitudinal pressure Pl{x). In 
terms of the electric and magnetic field strengths squared, 
as defined in Eq. ( 56 1 , they are given by 



e(x) 
P T (x) 



ce- 
f 2 



f 2 



P L { X )=C E ^- 



f 2 



CB 



cb 



B 2 T 



B 2 T 



(78) 
(79) 
(80) 



SEf(x) = f'{r,) e?(x x ) , (75) 
6E$(x) = 2ra" x /(„) [/'(t?)]" 1 £ , ^ 



1=1,2 



where /'(?7) is the lattice derivative of the function f(rj) 
according to 



f'(v)=a- 1 {f(r 1 )-f(ri-v)} 



(77) 



and the Gauss constraint is implemented explicitly in 
Eqns. f75| and d76l). 



These quantities are gauge invariant and satisfy the re- 
lation e = 2Pt + Pl at every position in space and time. 
In addition to the above quantities, we also study equal 
time correlation functions 

IL 2 L (T,v,v') = {PL(T,x^r 1 )P L {r,y ± , V ')} T , (81) 

where (.)t denotes ensemble average and average over 
transverse coordinates. We focus on the correlator in 
Fourier space with respect to relative rapidity, which is 
given by 



n£(T,i/) = L~ l / d v drf U l (t,v,v') e 



-iu(rj-rj') 



(82) 



and usually show results for Ili(^), i.e. the square root 
of the above expression. 
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